
#ifndef MATRIX_HPP
#define MATRIX_HPP

// A Matrix is a vector of vectors
//
// $Id: Matrix.hpp 2 2012-05-22 17:55:16Z gerry@gaboury.biz $

#include "Vector.hpp"

using namespace std;

template <class T> class Matrix : public Vector < Vector < T > >
{

private:

  // disable these methods
  virtual void pop_back () { }
  // row & col
  int m_row, m_col;
  ios::fmtflags m_ff;

protected:

  bool assertValid()
  {
    if ( this->size() != m_row )
      return false;
    for (int row=0; row < m_row; ++row )
      {
	if ( (*this)[row].size() != m_col )
	  return false;
      }
    return true;
  }
  //
  // update after using push_back
  //
  void update()
  {
    m_row = this->size();
    if ( m_row > 0 )
      {
	m_col = (*this)[0].size();
      }
  }

public:

  // access dimensions
  virtual int row() const { return m_row; }
  virtual int col() const { return m_col; }
  virtual void IOS_Format ( ios::fmtflags p_ff )
  {
    m_ff = p_ff;
  }

  virtual void push_back ( const Vector<T>& p_vec )
  {
    Vector< Vector<T> >::push_back(p_vec);
    this->update();
  }

  // constructors
  Matrix () : Vector < Vector < T > > (), m_row(0), m_col(0), m_ff(ios::fixed) { }

  Matrix (int p_row, int p_col) : 
    Vector< Vector< T > > (p_row), m_row(p_row), m_col(p_col), m_ff(ios::fixed)
  {
    for (int row=0; row < p_row; ++row )
      (*this)[row].resize(p_col);
  }

  Matrix ( int p_row, int p_col, T p_val ) : 
    Vector< Vector< T > > (p_row), m_row(p_row), m_col(p_col), m_ff(ios::fixed)
  {
    for (int row=0; row < p_row; ++row )
      (*this)[row] = Vector<T> ( p_col, p_val );
  }

  Matrix ( int p_row, int p_col, T* p_array ) : 
    Vector< Vector< T > > (p_row), m_row(p_row), m_col(p_col), m_ff(ios::fixed)
  {
    T* ptr = p_array;
    for (int row=0; row < p_row; ++row )
      {
	(*this)[row] = Vector<T> ( p_col );
	for ( int col=0; col < p_col; ++ col )
	  {
	    (*this)[row][col] = *p_array;
	    ++p_array;
	  }
      }
  }

  // copy constructor
  Matrix ( const Matrix<T>& p_mat ) :
    Vector < Vector < T > > ( p_mat.row() ),
    m_row ( p_mat.row() ), m_col ( p_mat.col() ), m_ff(ios::fixed)
  {
    for (unsigned row = 0; row < m_row; ++row )
      {
	(*this)[row] = p_mat[row]; // copy the row vectors
      }
  }

  // test for symmetry
  bool virtual isSymmetric() const
  {
    if ( m_row != m_col )
      return false;

    for (int row=0; row < m_row; ++row )
      {
	for ( int col=0; col < row; ++ col )
	  {
	    if ( (*this)[row][col] != (*this)[col][row] )
	      return false;
	  }
      }
    return true;
  }

  // empty matrix
  Matrix& zero()
  {
    for ( int i = 0; i < m_row; ++i )
      {
	for ( int j = 0; j < m_col; ++j )
	  (*this)[i][j] = 0;
      }
    return *this;
  }
  // identity matrix
  Matrix& identity()
  {
    this->zero();
    for ( int i = 0; i < m_row; ++i )
      {
	(*this)[i][i] = 1;
      }
    return *this;
  }
  // diagonal matrix
  Matrix& diagonal ( T p_val )
  {
    this->zero();
    for ( int i = 0; i < m_row; ++i )
      {
	(*this)[i][i] = p_val;
      }
    return *this;
  }
  // diagonal matrix
  Matrix& diagonal ( const Vector<T>& p_vec )
  {
    if ( p_vec.size() != m_row )
      throw (range_error("Matrix::diagonal(Vector) : dimension error"));
    this->zero();
    for ( int i = 0; i < m_row; ++i )
      {
	(*this)[i][i] = p_vec[i];
      }
    return *this;
  }

  // const Vector<T>& operator[] ( int row ) const
  // {
  //   if ( row >= m_row )
  //     throw (range_error("Matrix::operator[] : out of range."));
  //   return (*this)[row];
  // }

  Matrix<T> operator+ ( const Matrix<T>& p_mat ) const
  {
    if ( p_mat.row() != m_row || p_mat.col() != m_col )
      throw (range_error("Matrix::operator+ : Matrix dimension mismatch"));

    Matrix<T> Mtmp ( m_row , m_col );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < (*this)[row].size(); ++col )
	  {
	    Mtmp[row][col] = (*this)[row][col] + p_mat[row][col];
	  }
      }
    return Mtmp;
  }

  Matrix<T> operator- ( const Matrix<T>& p_mat ) const
  {
    if ( p_mat.row() != m_row || p_mat.col() != m_col )
      throw (range_error("Matrix::operator- : Matrix dimension mismatch"));

    Matrix<T> Mtmp ( m_row , m_col );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < (*this)[row].size(); ++col )
	  {
	    Mtmp[row][col] = (*this)[row][col] - p_mat[row][col];
	  }
      }
    return Mtmp;
  }

  Matrix<T> operator- ( const Vector<T>& p_vec ) const
  {
    if ( p_vec.size() != m_row )
      throw (range_error("Matrix::operator-(Vector) : size error."));

    Matrix<T> Mtmp ( m_row , m_col );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < (*this)[row].size(); ++col )
	  {
	    Mtmp[row][col] = (*this)[row][col] - p_vec[row];
	  }
      }
    return Mtmp;
  }

  Matrix<T>& operator+= ( const Matrix<T>& p_mat )
  {
    if ( p_mat.row() != m_row || p_mat.col() != m_col )
      throw (range_error("Matrix::operator+= : Matrix dimension mismatch"));

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < (*this)[row].size(); ++col )
	  {
	    (*this)[row][col] += p_mat[row][col];
	  }
      }
    return *this;
  }

  Matrix<T>& operator-= ( const Matrix<T>& p_mat )
  {
    if ( p_mat.row() != m_row || p_mat.col() != m_col )
      throw (range_error("Matrix::operator-= : Matrix dimension mismatch"));

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < (*this)[row].size(); ++col )
	  {
	    (*this)[row][col] -= p_mat[row][col];
	  }
      }
    return *this;
  }

  // matrix multiplication
  Matrix<T> operator* ( const Matrix<T>& p_mat ) const
  {
    if ( p_mat.row() != m_col )
      throw (range_error("Matrix::operator* : Matrix dimension mismatch."));

    // result
    Matrix<T> Mtmp ( m_row , p_mat.col() );

    for (unsigned irow = 0; irow < m_row; ++irow )
      {
	for (unsigned jcol = 0; jcol < p_mat.col(); ++jcol )
	  {
	    T tmp(0);
	    for ( unsigned k = 0; k < m_col; ++k )
	      tmp += (*this)[irow][k] * p_mat[k][jcol];

	    Mtmp[irow][jcol] = tmp;
	  }
      }
    return Mtmp;
  }

  // Kronecker Product
  Matrix<T> KronProd ( const Matrix<T>& p_mat ) const
  {
    // result
    Matrix<T> Mtmp ( m_row * p_mat.m_row , m_col * p_mat.m_col );

    for (unsigned lrow = 0; lrow < m_row; ++lrow )
      {
	for (unsigned rrow = 0; rrow < p_mat.m_row; ++rrow )
	  {
	    for (unsigned lcol = 0; lcol < m_col; ++lcol )
	      {
		for ( unsigned rcol = 0; rcol < p_mat.m_col; ++rcol )
		  Mtmp[lrow*p_mat.m_row+rrow][lcol*p_mat.m_col+rcol] = 
		    (*this)[lrow][lcol] * p_mat[rrow][rcol];
	      }
	  }
      }
    return Mtmp;
  }

  // scalar assignment
  const Matrix<T>& operator=( const T& p_scalar )
  {
    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    (*this)[row][col] = p_scalar;
	  }
      }
    return *this;
  }
  // scalar multiplication
  Matrix<T> operator*( const T& p_scalar ) const
  {
    // result
    Matrix<T> Mtmp ( m_row , m_col );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    Mtmp[row][col] = (*this)[row][col] * p_scalar;
	  }
      }
    return Mtmp;
  }
  // scalar division
  Matrix<T> operator/( const T& p_scalar ) const
  {
    // result
    Matrix<T> Mtmp ( m_row , m_col );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    Mtmp[row][col] = (*this)[row][col] / p_scalar;
	  }
      }
    return Mtmp;
  }

  // transpose
  //Matrix<T> operator!() const
  //Matrix<T> operator^(int i) const
  Matrix<T> Tr() const
  {
    // result
    Matrix<T> Mtmp ( m_col , m_row );

    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    Mtmp[col][row] = (*this)[row][col];
	  }
      }
    return Mtmp;
  }

  // scalar multiplication
  Matrix<T>& operator*=( const T& p_scalar )
  {
    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    (*this)[row][col] *= p_scalar;
	  }
      }
    return *this;
  }
  // scalar division
  Matrix<T>& operator/=( const T& p_scalar )
  {
    for (unsigned row = 0; row < m_row; ++row )
      {
	for (unsigned col = 0; col < m_col; ++col )
	  {
	    (*this)[row][col] /= p_scalar;
	  }
      }
    return *this;
  }

  // vector multiplication
  Vector<T> operator* ( const Vector<T>& p_vec ) const
  {
    if ( p_vec.size() != m_col )
      throw (range_error("Matrix::operator*(Vector) : dimension error"));

    // result
    //VectorStore* VS = VectorStore::instance();
    //Vector<T>* Vtmp = VS->newVector<T>( m_row );
    Vector<T> Vtmp ( m_row );

    for (unsigned row = 0; row < m_row; ++row )
      {
	T tmp(0);
	for ( unsigned k = 0; k < p_vec.size(); ++k )
	  tmp += (*this)[row][k] * p_vec[k];

	Vtmp[row] = tmp;
      }

    return Vtmp;
  }

  void write (ostream *os, int p_prec=5, int p_width=10 ) const
  {
    // set output format
    ios::fmtflags oldopt = os->setf( m_ff );
    os->precision(p_prec);

    for (int i=0; i < this->size(); ++i)
      (*this)[i].write(os,p_prec,p_width);
    *os << endl;

    // reset output format
    os->flags( oldopt );
  }

  void read (istream *is)
  {
    // Vector read
    for (int i=0; i < this->size(); ++i)
      (*this)[i].read(is);
  }
};

template <class T> ostream &operator<< ( ostream &os,Matrix<T> &p_matrix )
{
  p_matrix.write(&os);
  return os;
}
template <class T> ostream &operator<< ( ostream &os,const Matrix<T> &p_matrix )
{
  p_matrix.write(&os);
  return os;
}
template <class T> istream &operator>> ( istream &is,Matrix<T> &p_matrix )
{
  p_matrix.read(&is);
  return is;
}

#endif // MATRIX_HPP
